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ESTIMATION OF EFFECTIVE HYDROLOGIC PROPERTIES OF SOILS 


FROM OBSERVATIONS OF VEGETATION DENSITY 


by 


tol:n e. tellers 

and 

PETER S. EAGLESON 


ABSTRACT 


An existing one-dimensional model of the annual water balance 
is reviewed. Slight improvements are made in the method of calculating 
the bare soil component of evaporation, and in the way surface retention 
is handled. A natural selection hypothesis, which specifies the equilib- 
rium vegetation density for a given, water-limited, climate-soil system, 
is verified through comparisons with observed data and is employed in the 
annual water balance of watersheds in Clinton, Ma., and Santa Paula, Ca., 
to estimate effective areal average soil properties. Comparison of CDF’s 
of annual basin yield derived using these soil properties with observed 
CDF’s provides excellenc verification of the soil-selection procedure. 
This method of parameterization of the land surface should be useful with 
present global circulation models, enabling them to account for both the 
non-linearity in the relationship between soil moisture flux and soil 
moisture concentration, and the variability of soil properties from place 
to place over the earth’s surface. 


2 


ACKNOWLEDGEMENTS 


This work was completed with the support of the National 
Aeronautics and Space Administration (NASA) under Grant NSG 5306. 

The work was performed by Tobin Tellers, research assistant in 
Civil Engineering at MIT, and constitutes his thesis presented in partial 
fulfillment of the requirements for the degree of Master of Science in 
Civil Engineering. This study was supervised by Peter S. Eagleson, 

Professor of Civil Engineering, who also provided the theoretical background 
material on which it is based. 

Thanks also are due to Bernhard Metzger, graduate student, who 
provided many helpful suggestions and comments, and to Anne Clee, who 
performed all the necessary typing. 


TABLE OF CONTENTS 


ABSTRACT 2 

ACKNOWLEDGEMENTS 3 

TABLE OF CONTENTS 5 

LIST OF PRINCIPAL SYMBOLS 6 

LIST OF TABLES 11 

LIST OF FIGURES 12 

Chapter 1 INTRODUCTION 14 

Chapter 2 OBJECTIVES 15 

Chapter 3 REVIEW OF LITERATURE 16 

Chapter 4 THEORETICAL BACKGROUND 23 

4.1 The Water Balance 23 

4.2 Evapotranspiration 26 

4.3 Surface Runoff 36 

4.4 Vegetal Equilibrium Hypothesis 40 

Chapter 5 METHOD OF APPROACH 53 

5.1 Vegetal Equilibrium Hypothesis 53 

5.2 Estimation of Effective Soil Properties 59 

Chapter 6 PRESENTATION OF RESULTS 69 

6.1 Verification of Vegetal Equilibrium 

Hypothesis 69 

6.2 Estimation of Effective Average Areal 

Soil Properties 71 

Chapter 7 SUMMARY AND CONCLUSIONS 91 

5 

m *Jt 




Chapter 8 RECOMMENDATIONS FOR FUTURE WORK 94 
REFERENCES 95 
Appendix A 98 
Appendix B 110 


6 




LIST OF PRINCIPAL SYMBOLS 


A 

A 


o 


c 


d 




e 

* 


f . 
1 


G 

h 

h 


albedo 

gravitational infiltration rate as modified by capillary 
rise from water table, cm/sec. 

pore disconnectedness index 

diffusivity iftdax 

long-term average rate of potential evaporation, cm/sec. 
long-term average rate of actual evaporation, cm/sec, 
rate of transpiration, cm/day 
exfiltration parameter 

average annual potential evapotranspiration, cm 
annual surface retention, cm 

soil moisture evaporation from bare soil fraction of 
surface, cm 

total evapotranspiration, cm 

annual evapotranspiration, cm 

annual evaporation from soil moisture, cm 

transpiration from vegetated fraction of surface, cm 

exfiltration rate, cm/day 

exfiltration capacity, cm/day 

infiltration rate, cm/day 

infiltration capacity, cm/day 

gravitational infiltration parameter 

storm depth, cm 

surface retention capacity, cm 


7 


H 


i 


k(l) 

K(l) 


k 


v 


L 

e 

in 


"h 

m i 

% 

\ 

m 

t 

i 

m 

V 

m 

T 

M 

M 

o 

n 

N 



average residual sensible heat flux, ly/min 

precipitation rate, cm/day 

saturated intrinsic permeability, sq. cm 

saturated effective hydraulic conductivity, cm/sec. 

plant coefficient 

latent heat of vaporization 

pore size distribution index 

mean storm depth, cm 

mean rainfall intensity, cra/day 

average annual precipitation, cm 

mean time between storms, days 

mean storm duration, days 

mean number of storms per year 

mean length of rainy season 

vegetal canopy density 

equilibrium vegetal canopy density 

medium effective porosity 

percent cloud cover 

annual precipitation, cm 

average rate of net outgoing long wave radiation, ly/min. 

average rate of insolation, ly/min. 

annual groundwater runoff, cm 

annual surface runoff, cm 

annual x*ainfall excess, cm 

time and spatial average soil moisture concentration 
in surface boundary layer 


8 


s 

s 



w 


a 

$ 


w 


n 


K 

x 


w 


w 




relative humidity 

1/2 

exfiltration sorptivity, cm/s 
time between storms, days 

duration of exfiltration rate at exfiltration capacity, days 

duration of exfiltration at the rate ip, days, 
time at which surface retention reaches saturation 
during precipitation, days 

storm duration, days 

normal annual temperature, °C 

upward apparent pore fluid velocity representing capillary 
rise from the water table, cm/sec. 

annual yield, cm 

reciprocal of average rainstorm intensity, days/cm 

reciprocal of average time between storms, days ^ 

3 

specific weight of water, dynes/cm 
reciprocal of average storm duration, days * 
reciprocal of mean storm depth, cm * 
parameter of Gamma distribution of storm depth 
parameter of Gamma distribution of storm depth, cm ^ 
dynamic viscosity of water, poises 

3 

mass density of evaporating water, gr/cm 
capillary infiltration parameter 
surface tension of water, dynes/cm 
latitude 

dimensionless exfiltration diffusivity 
dimensionless infiltration diffusivity 


9 


$(m) 

pore shape parameter 

H'(l) 

saturated soil matrix potential 

y/a 

atmospheric parameter 

E[ ] 

expected value of [ ] 

J( ) 

evapotranspiration function 

r( ) 

Gamma function 

Y(a, x) 

incomplete Gamma function 


(suction) 


10 


LIST OF TABLES 


Table 

Title 

Page No 

4.1 

Independent Parameters of Representative 
Soils 

39 

5.1 

Albedo of Natural Surfaces 

57 

6.1 

Input Climate and Vegetation Parameters 

73 

6.2 

Annual Water Balance Components, Santa Paula, 
Ca. Estimated Soil Properties. Equation 
(5.6), M q - .4 

74 

6.3 

Annual Water Balance Components, Santa 
Paula, Ca. Silty-loam Soil Properties. 
Equation (5.6), M q » .4 

75 

6.4 

Annual Water Balance Components, Clinton, 
Ma. Estimated Soil Properties. Equation 
(5.6), M q - .8 

77 

6.5 

Annual Water Balance Components, Clinton, 
Ma. Silty-loam Soil Properties. Equation 
(5.6), M q * .8 

78 

6.6 

Annual Water Balance Components, Santa 
Paula, Ca. Estimated Soil Properties. 
Equation (5.5), M q ■ .4 

81 

6.7 

Annual Water Balance Components, Clinton, 
Ma. Estimated Soil Properties. Equation 
(5.5), M q - .8 

83 

6.8 

Annual Water Balance Components, Santa 
Paula, Ca. Estimated Soil Properties. 
Equation (5.6), M q * .2 

85 

6.9 

Annual Water Balance Components, Clinton, 
Ma. Estimated Soil Properties. Equation 
(5.6), M q = .6 

87 

6.10 

Annual Water Balance Components, Clinton, 
Ma. Estimated Soil Properties. Equation 
(5.6), M q = .9 

88 


11 


LIST OF FIGURES 


Figure 

4.1 

4.2 

4.3 

4.4 

4.5 

4.6 

4.7 

4.8 

4.9 

4.10 

4.11 

5.1 

5.2 

5.3 

5.4 


Title 


Page No. 


Schematic Representation of Vegetated Soil 
Column during an Interstorm Period 

Interstorm Evaporation from Bare Soil 


27 

30 


Integration Regions for Calculation of 
Expected Value of Interstorm Evapotranspiration 

Surface Runoff Generation during Typical Storm 

(t > t ) 
r o 

Sensitivity of Mean Annual Soil Moisture to 

Vegetal Canopy_Density for Typical Soils 

(k =1 and w/e « 1) 
v p 

Sensitivity of Mean Annual Evapotranspiration 

to Vegetal Canopy Density for Typical Soils 

(k = 1 and w/e « 1) 
v p 

Evapotranspiration Function for Natural Systems 
and Equilibrium Vegetal Canopy Density (k = 1, 
= 0.0, K = .5, w/e^ « 1) X 

Evapotranspiration Function for Natural Systems 

(k = 1, h = 0.0, k * .5, w/e << 1) 
v o p 

Sensitivity of Evapotranspiration Function to 
8h / e and Ah (k * 1, K = .5, w/e « 1) 

o p O V p 

Sensitivity of Evapotranspiration Function to K 

(8h /e = .lj k = 1, w/e « 1) 
o p J v ’ p 

Sensitivity of Evapotranspiration Function to 

Bh /e and k (k * .5, Ah * .05, w/e << 1) 
o p V * o p 

Insolation at Earth’s Surface 

Long-Term Average Potential Evapotranspiration 


31 

37 

41 

43 

46 

47 

49 

50 

51 
56 
58 


Saturated Permeability vs. Pore Size Distribution 
Index 


60 


Water Balance Solutions Using Soil Properties 
from Equation (5.3) 


67 


12 


Figure 

Title 

Page Ho. 

6.1 

Verification of Vegetal Equilibrium Hypothesis 

70 

6.2 

Frequency of Annual Basin Yield, Santa Paula, Ca. 

76 

6.3 

Frequency of Annual Basin Yield, Clinton, Ma. 

79 

6.4 

Sensitivity of Annual Basin Yield to Two 
Methods of Handling Surface Retention, Santa 
Paula, Ca. 

82 

6.5 

Sensitivity of Annual Basin Yield to Vegetal 
Canopy Density, Santa Paula, Ca. 

86 

6.6 

Sensitivity of Annual Basin Yield to Vegetal 
Canopy Density, Clinton, Ma. 

90 


Chapter 1 


INTRODUCTION 


In order to increase the accuracy of global climate models, 
a more sophisticated representation of the land surface boundary condi- 
tion is required than that which is presently employed (GARP, 1978). The 
interaction, in particular the water flux, between the atmosphere and 
the soi v -vegetation system at this boundary is highly non-linear in 
nature, and is not simply defined. Any attempt to satisfactorily 
account for this non-linearity in a model must incorporate two effects 
which are not included in current models: 

1. variability of soil properties and soil moisture 
dynamics from place to place over the earth’s surface, and 

2. non-linearity in the relationship between soil moisture 
flux and soil moisture concentration. 

In this work, it is intended to make use of a one-dimensional 
water balance at the land-air interface in order to parameterize the 
climate-soil-vegetation relationship in such a way as to reflect the 
non-linearity and areal variability. 
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Chapter 2 


OBJECTIVES 

The specific objectives of this work are twofold: 

The first objective is verification of a vegetal equilibrium 
hypothesis developed by Eagleson (1978f). This hypothesis proposes 
that the natural vegetation density it; a watershed will seek, through 
natural selection, an optimal "climax" value at which available soil 
moisture is a maximum. Comparison of a theoretical curve of evapo- 
transpiratlon versus canopy density based on this hypothesis with 
observed data will provide the necessary check on the accuracy of the 
hypothesis. 

The second objective is establishment of an algorithm for 
estimating effective areal soil properties from observations of 
vegetation density by using the natural selection hypothesis in a 
one-dimensional water balance model. By defining the level of evapo- 
transpiration from soil moisture through observations of the canopy 
cover density, it may be possible, knowing the climate, to determine 
the soil properties that enable the soil-vegetation system to respond 
at the Indicated level. The estimated values of these parameters can 
then be used in the water balance equation to evaluate desired compo- 
nents of the water flux. Verification of the desired algorithm will 
be sought through comparison of computed and observed statistics 
of annual yield. 
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Chapter 3 


REVIEW OF LITERATURE 

Past efforts to model the coupling among physical processes 
of the atmosphere, soil, and vegetation across the land-surface inter- 
face have been largely of two types: 

1. Numerical studies which employ detailed formulations 
of the processes involved. Examples of such studies are those of 
Philip (1957), Sasamori (1970), Deardorff (1977), and Philip and de Vries 
(1957) . Although these models simulate the system response to climatic 
inputs very well, they usually do so in terms of a large number of 
climate, soil, and vegetal parameters. Due to their complexity and 

the detailed data requirements for their validation, these studies 
have found little application in general circulation models. 

2. Empirical studies which utilize validated interrelation- 
ships among the principal variables. Because of the ease of their 
application, and negligible programming and data requirements, most 
global climate models use this type of parameterization of the land- 
surface boundary with regard to actual evapotranspiration, average 
soil moisture content, and runoff. 

The primary GCM's today utilize the approach first introduced 
by Manabe (1969) to parameterize the land surface boundary condition. 

In this approach, the above mentioned parameters are handled in the 
following way: 
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A. Evapotransplratlon 

Actual evapotransplratlon is related to potential evapo- 
transplration linearly through the soil moisture and a single soil 
parameter following the work of Budyko (1956) . This parameterization 
is 


r 

e T s/k, 8 £ k 

~ » * (3.1) 

P 1 , k < s < 1 

in which 

■ actual rate of evapotransplratlon 

e p ■ potential rate of evapotransplratlon 

s - effective soil moisture concentration 

k ■ empirical coefficient, 0<k£ 1 generally assumed 

to be constant everywhere 


As mentioned above, the only soil parameter appearing in this 
model is the empirical coefficient, k. This representation grossly 
distorts the sensitivity of e T to s and makes no allowance for the spatial 
variance of this sensitivity due both to soil type and to the presence 
of vegetation. 

More recently, Lettau (1969) and Lettau and Baradas (1973), 
in their "evaporation climatonomy" formulation, refine the water 
balance evapotransplratlon term through use of an energy balance. This 
approach seeks theoretical solutions in the form of "response functions" 
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(i.e. evapotranspiration cycles, temperature, etc.) as a physical con- 
sequence of a mathematically defined "forcing function" of the environ- 
mental system. However, parameterization is achieved without any 
explicit consideration of the soil and vegetal properties which will 
control the evaporation under all but the most humid conditions. The 
only input parameters linked to the land surface are evaporivity, e* 
which is a non-dimensional measure of the capacity of land surfaces to 
utilize solar energy for the evaporation of rainfall received in a 
specified time interval, and t* which denotes a characteristic soil 
moisture residence time. Values for these parameters are either 
assumed on the basis of empirical data, or are estimated from a systematic 
classification of watersheds according to morphology, soil structure and 
permeability, vegetation cover, etc. The lumping of all these para- 
meters into a single term in no way fully represents the complex 
interrelationships between the various processes involved in the 
water balance. 

Other studies concerning the evapotranspiration term are 
those of Czarnowski (1964) and Ritchie (1972), and Ritchie and Burnett 
(1971) . Czarnowski assumes that total evapotranspiration is a sum of 
plant transpiration, evaporation from surface retention, and evaporation 
from soil, and that these values are functions of vegetation density, 
and consequently of climatic factors. He treats the development of 
plant cover as a function of the form 

_P_ 

V 

M = 1 - e “ (3.2) 
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where 

P « precipitation, mm 

V m * sum of mean daily deficits of air humidity, ram Hg 
Finally, he concludes that evapotranspiration can be expressed as 


V - V M 
m 


1.17M + 


.4 

/M - 


or 


M 


m 


1.17M + 


♦ 4 

/F 


(3.3) 


(3.4) 


where the constants, 1.17 and .4 are determined by a least squares fit 
to empirical data obtained primarily from cultivated agricultural 
lands . 

Ritchie (1972) and Ritchie and Burnett (1971) develop a set 
of empirical functions relating leaf area index and fractional net 
radiation at a soil surface for a row crop to plant evaporation 
efficiency. These equations may be written 


„ t> 398L , ,, 

R = R e Ai (3.5) 

ns no 

— * -70L 1/2 - .21 (3.6) 

e p Ai 

where 

L a , - leaf area index 
Ai 


19 


net radiation reaching soil surface 


R no ~ net rac *i ation above plant canopy 

Relating leaf area index to canopy density using Equation (3.5) and 
the assumption that 


R 

= 1 - M (3.7) 

no 


gives 


J 



ln(l-M) 


1/2 



.21 


(3.8) 


Again, the constants appearing in the above equations are determined 
from the method of least squares. 

While both of the above formulations are attempting to relate 
evapotranspiration to more physically-significant parameters, there is 
little inclusion of the actual physics. Since a linear regression is 
performed to obtain the above equations, there is a lack of generality 
and understanding of the sensitivity to other parameters besides vege- 
tation 


B. Soil Moisture 

The change in the average soil moisture concentration is 
determined from a water balance relation written 
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(3.9) 


U 3s 
nh St 


1 ■ e T " Y S " Y Ci 


in which 

h * thickness of surface layer 

i * intensity of rainfall 

Yg ■ intensity of surface runoff 

Yg * intensity of percolation of water out of the surface 
layer to groundwater 

The product, nh, represents the maximum water content of 
the surface layer and is assigned a value which is common for all 
soil surfaces. 

C. Runoff 

Runoff, as written in Equation (3.9) consists of two different 
components. Surface runoff is regulated by the infiltration of rainfall 
and additions to soil moisture. Groundwater runoff is governed by the 
state of soil moisture concentration. All global climate models use 
highly simplified empirical formulae which lump these two dynamically 
different runoff-generating processes into total yield relations of 
the form 


Y = Y(i,e T , s) 
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These relations include one or more coefficients which may incorporate 
spatial variability, but there is no physical basis for their selection 
without natural yield measurements. 

The models referred to in the preceding paragraphs include 
those of Arkawa (1972, U.C.L.A.); Somerville et al. (1979, G.I.S.S.); 
Gates and Schlesinger (1977, Rand-O. S. U. ) ; Sellers (1973, Arizona); 
and Corby et al. (1978, B.M.O.). In all of these models, there is 
no use of the present high level of physical understanding of the 
natural processes involved to develop a generalized, accurate repre- 
sentation of the land-surface interface. 

Eagleson (1978a,b,c,d,e,f ,g) , has developed a generalized 
water balance based upon simplified physics of the component processes. 
The development is sufficiently rigorous to capture the essential 
system dynamics yet simple enough to permit analytical solution. 

The model produces valuable insights into the interactive role of 
soil moisture in the determination of climate. Foremost in this 
development is the accounting for the areal variability of soil pro- 
perties over the earth’s surface and the reflection of the inherent 
non-linearity in the relationships between soil moisture concentration 
and the interfacial moisture fluxes. This model will be presented 
and utilized in the following chapters to attain the objectives 
stated in Chapter 2. 
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Chapter 4 


THEORETICAL BACKGROUND 


4. 1 The Water Balance 

The major source of theoretical background used here is the 
work of Eagleson (1978a, b, c, d, e, f ,g) . In these papers, a one- 
dimensional water balance based on soil moisture dynamics and 
statistics of climatic data is derived. This water balance, expressed 
in terms of annual expected values, may be represented as 

E[E p ] - E[E X ] + E[R g ] + E[R G j (4.1) 

A AAA 

and 

- E[R g ] + E[R g ] (4.2) 

A A 


■ expected value of I ] 

* annual precipitation 

* annual evapotranspiration 

■ annual surface runoff 

* annual groundwater runoff 

■ annual yield 


where 


E[Y k ] 


E[ ] 
P. 
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An analytic expression is obtained for Equation (4.1) by deriving the 
individual components through the use of derived probability distribu- 
tions and one-dimensional dynamic equations approximating the physics 
of the separate soil moisture fluxes. These expressions are then 
introduced into Equation (4.1) to produce the equation for the (soil 
moisture) water balance 

m_ (l-e' G_2a r (a + 1 )a'°) - 

A 

E[E p ] J(E,M,k v ,h o ) - E[E r ] + m t K(l) s o C - TW 


for 


E[E )/mp <e G2£7 r(cr + l)a 0 (4.3a) 

A A 

(The term to the left of the equal sign is infiltration, the first 
term to the right, total evapotranspiration, the second is evaporation 
from surface retention, and the last two terms are groundwater runoff 
[the first is groundwater recharge and the last is groundwater loss]). 
Otherwise, 


mp = E[E p ] J(E,M,k v ,h o ) + m T K(l)s o C - TW (4.3b) 

A A 

In the above, 

E = annual surface retention 
r A 
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average annual potential evapotranspiration 



J 

G 


a 

E 

M 

k 

v 


h 


o 


s 

o 


K(l) 


T 


w 


% 


evapotranspiration efficiency 
gravitational infiltration parameter 
capillary infiltration parameter 
evapotranspiration parameter 
vegetation canopy density 
plant transpiration coefficient 
mean length of rainy season 
surface retention capacity 
average annual soil moisture 
saturated hydraulic conductivity 
1 year, seconds 

apparent velocity of capillary rise 
mean annual precipitation 
pore disconnectedness index 


It will be helpful and important to review the development of 
the expressions for evapotranspiration and surface runoff, and to 
present an alternative approach for the former, and a slightly different 
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interpretation for the latter. 


4.2 Evapotranspiration 

The expected value of annual evapotranspiration is derived 
(Eagleson, 1978d) by calculating bare soil evaporation and vegetal 
transpiration for an interstorm period as functions of properties of 
the storm sequence, the surface, the soil, and the average rate of 
potential evapotranspiration, using observed distributions of the 
random climatic variables , and averaging over the rainy season. The 
bare soil evaporation and plant transpiration are determined by con- 
sidering the vertical flux of moisture in a soil column. In Figure 
4.1, the modeled column of soil and the different moisture fluxes 
are sketched. In this figure 

f = bare soil exfiltration rate 
e 

M = vegetation canopy density 

e v = vegetation transpiration rate 

K(0 q ) = effective hydraulic conductivity at long-term average 
soil moisture 

It is assumed here that 

1. Soil moisture throughout the surface boundary layer is 
spatially uniform at the start of each interstorm period 
at the long-term average value, s = s q ; 

2. The medium is semi-inf inite; 
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3. 


The vegetation is distributed uniformly, and its roots 
extend uniformly into the entire volume of the soil in 
the surface boundary layer. This implicitly assumes that 
the plant species have adapted by natural selection to 
a density and root structure which is in balance with 
the available soil moisture; 

4. The rate of moisture extraction by the roots is in 

equilibrium with the transpiration rate by the leaves, 
thus forming a uniformly distributed sink for soil 
moisture of strength, Me v » 


Following the work of Philip (1969), Eagleson writes the 
total decrease in soil moisture during infiltration: 


(0 - 0)dz 

o 

0 


zd0 = F (t) + [K(0 ) + Me ] t 
. e o v 

e i 


(4.4) 


where F @ (t) is the cumulative exfiltration in centimeters. 

The integral on the left-hand side is evaluated in the manner 
of Philip (1960) . Assuming a vertical flow passage of constant 
cross-section, the exfiltration rate is found to be 

_1 

f (t) = T St 2 ~ " Me v (4 * 5) 

e Z e z ± o v 


Note that this neglects the restriction, by vegetation canopy density, 
of the bare soil area through which exfiltration occurs. Further 
simplification and analysis result in the exfiltration capacity: 
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(A. 6) 


Vi 


* 


f 

e 




-Me 

v 


where S is the exfiltration sorptivity. 

© 

A typical interstorm period, and the relationship between the 

various rates of evaporation and time for bare soil is illustrated 

in Figure 4.2. In this figure, e p is the potential rate of bare 

soil evaporation, which is considered a constant. The times, t Q and 

t , are evaluated by assuming 
e 

e « f* (4,7a) 

p e 


when 


£*(Ve) 


feVo* 


(4.7b) 


and that 


1* « 0 (4.8a) 

e 


when 

t - t (4.8b) 

e 

respective] v. Exfiltration capacity, and the times t fi and t Q , are 
then used by Eagleson along with the relationships represented in 
Figure 4.2 to determine total evapotranspiration, E T * To do this. 
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E x from a unit land surface is proportioned according to 

E t - (l-M)E g + ME y (4.9) 

in which 

E s = bare soil evaporation from soil moisture plus eva- 
poration from soil surface retention 
E v * evapo transpiration from vegetation plus evaporation 
from plant surface retention. 


It is not necessary to present the development of E[E ] here. 

s 

This is done by calculating the volume under the solid line in Figure 
4.2, multiplying by the joint probability distribution of storm depth 
and time between storms, and integrating over the regions shown in 
Figure 4.3. What is important to note is the previously mentioned 
approximation made in the development of the bare soil exfiltration 
capacity. The expression obtained for E[E ], bare soil evapotrans- 

S j 

piration for one interstorm period, from the above procedure is 


E[E ] = 
3 


e 

_R 

3 


Y[K,Ah Q ] 

t(k) 


-K 


3h /tt " 

1+ -°— E 


Ah 


O -J 


Yt<,Ah +3h /e ] -BE 

*— e 

TOO 


1- 


y[K,xh o ] 

t(k) 


-BE-3h fe 

o p 


1-e 


[14Mk v + (2B) 1/2 E-w/e p ] 
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-CE-Bh /e , /9 

+ e P [Mk y + (2C) ' E - w/oJ + 


(2E) 


1/2 e B q/ P [Y(f,CE) - Y(f»BE)| 


Bh/e 

i + — 2 — E 

1 Ah 


-t — K 


Y [ K ,Ah o +3h o /V 

V(k) 


|(2E) 1/2 [ Y (|,CE)- Y (|,BE)] 


+ e 


-CE 


[Mk v + (2C) 1/2 E - w/e p ] - e _BE [Mk y + (2B) l/2 E - w/e p ] J 


(4.10) 


Here 


B 


1-M _ 
l*k v -w/e p 


+ 


M 2 k +(l-M)w/e 
v p 

2(1+Mk -w/e ) 2 
v p' 


(4.11) 


and 


C = T (Mk v ' w/e p ) 


(4.12) 


Also, 

B = reciprocal of mean time between storms 
A = parameter of Gamma distribution of storm depth 
k = paramater of Gamma distribution of storm depth 
h Q = surface retention capacity 
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Upon studying Figure 4.1 and Equations (4.4) and (4.6) 

it can be seen that the term, £ , is defined as the exfiltration rate 

e 

for bare soil, and F as the total volume of moisture exfiltrated from 

e 

the soil column across the bare soil surface. The rate, f , is 

e 

obtained by differentiating the volume, F e » with respect to time. 

The result of the differentiation leaves f multiplied by the area 

of bare soil. Thus, in the two-dimensional problem which includes 

the presence of vegetation, f should be multiplied by the term, 1-M, 

e 

to account for the fact that only a fraction of the land surface, the 
bare soil fraction, is exfiltrating at this rate. Equation (4.6) 
should be rewritten as 

(l-M)f* = t" 1/2 - Me (4.13) 

e 2 e v 


The new form for the expected value of bare soil evaporation, E , may 

s 

then be evaluated in terms of this altered expression for exfiltration. 


The new expression for E[E ] is 


’j 


E[E s j I ■ ~£ | r<K) 


vtMig - B e 

- e 


1 + 


B'n fe 

0 -P 
x h 


Y[K,Xh Q +Bh o /e p ] 

F(tc) 


1 - 


y[»c,Xh o ] 

TOO” 


1-e 


-BE-gh /e 
o p 


1+ 


[' 


(2B) 1/2 E - — +Mk 

■S’ V 

P 


-e 


-CE-3h /e 
o p 


1-M 


- Mk - (2C) 1/2 E 
^ v 


P 


1-M 


, /0 ~gh /e 
(2E) 1/2 a ° p 


[y(|,ce) - y(|,be)1 
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Sh/e 

1 +— 2— E. 


Xh 


— K 


Y[K, Xh o + ghp/%] 

r(K) 


-BE 


1 fw_ 
1-M U " 


Mk - (2B) 1/2 E 
v 


-CE| 1 f w 


Mk - (2C) 1/2 E 
v 


1-M 


(2E) 


V2 Ul 


2 , CE)-Y(f, BE) 


(A. 14) 


where 


B » 1/(1+Mk v ) 

C + 1/2 (Mk -w/e ) -2 
v p 


(4.15) 

(4.16) 


E j is obtained in the same way as before; by multiplying the bare soil 
term by (1-M), and the vegetation term by M. The result of this altera- 
tion on the expected value of annual basin evapo transpiration will be 
presented in a later section. Although this approach may seem more 
accurate than the original, its use will create other, and possibly 
greater problems. Attempting to expand the problem into two dimensions 
at this point will cause some inconsistencies concerning evaluation of 


the Philip exfiltration equation. This equation is 
0 

o 


zde - Ajt 1 / 2 + A 2 t + A 3 t 3/2 + . . . 


(4.17) 


e. 


Since this was developed for a one-dimensional formulation, the 
expressions obtained for the constants, A^, on the right hand side will 
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not necessarily apply to the two-dimensional situation. This can be 
seen by noting that in Equation (4.6), as M approaches the value of 1, 
the right hand side does not go to zero, as it should for a fully 
vegetated surface where there would be no bare soil exfiltration of 
soil moisture. So, although there are certain misgivings about 
Eagleson’s original derivation, the alternate approach presented above 
may involve more serious inaccuracies. However, for areas with a large 
vegetal canopy density, where the effect of the vegetation on bare soil 
exfiltration is large, this approach may come closer to reality than 
the previous one. 

4.3 Surface Runoff 

To derive the probability of storm surface runoff, Eagleson 
(1978) integrates the difference between rainfall intensity and the 
Philip infiltration equation over the duration of a rainstorm. Infil- 
tration is assumed to occur uniformly over both bare soil and vegetated 
portions of the surface. Illustrated in Figure 4.4 is a sequence of 
surface states beginning from t = 0 at the start of the rainfall period. 
In this figure 

h * surface retention capacity 

0 

1 = rainfall intensity 
* 

= infiltration capacity 

t = storm duration 
r 

A q = gravitational infiltration rate as modified by capillary 
rise from the water table 
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Initially, there is a withdrawal of rainfall to satisfy surface retention. 
If t is greater than h Q /i, as shown, this surface retention reaches its 
capacity, h Q . If t £ h Q /i, there would be no infiltration or runoff, 
and the surface retention would equal the storm depth, h. For the case 
illustrated in Figure 4.4, however, infiltration will begin at time 

•fa 

t = h /i. From this time until t » h /i + t , when f. * i, infiltration 
o o o i 

will take place at the rainfall rate, i. After this time, the capacity 

of the soil to infiltrate moisture is no longer larger than the rainfall 

intensity, and the excess is represented by the shaded area of the 

* 

figure. Rainfall excess, R , is then generated until time t = t . The 

S j r 

expected value of the rainfall excess is obtained in a manner similar to 

that of the evapo transpiration. 

A question may be raised relating to the hackling of the 

surface retention. In his development, Eagleson argues that the surface 

retention must be subtracted from the rainfall excess, since it is 

moisture that is not infiltrated into the soil. The expression he 

obtained for the expected value of annual rainfall excess is then 

E[R* ] = hl, [e“ G " 2a T(cr + l)a“ a ] (4.18) 

S A r A 

in which 

* 

R = annual rainfall excess 
S A 

The expected value of annual surface retention, E , is then 

r A 

subtracted from this to get the annual surface runoff. This charges the 
entire annual surface retention against those events producing rainfall 
excess, however. 
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CLAY 

CLAY-LOAM 

SILTY-LOAM 

k(D 

1.0xl0~ 10 

2.8xl0" 10 

1.2xlO -9 

n 

.45 

.35 

.35 

c 

12 

10 

6 


SANDY-LOAM 

2.5xl0~ 9 


.25 

4 


Table 4.1 

REPRESENTATIVE SOIL PROPERTIES 


39 


A slightly different interpretation of this rainfall excess 
results from a closer examination of Figure 4.4 and the shaded area 
therein. It is known, and mentioned above, that surface retention 
must be satisfied before any infiltration can occur. From this, it 
seems necessary to subtract the surface retention from the beginning 
of the rainfall period, as indicated in the figure, rather than from 
the rainfall excess at the end. The volume represented by the shaded 
area would then be equal to the surface runoff, and not surface runoff 
plus surface retention. The resulting water balance equation then 
becomes 

nip (1 - e" G_ 2 a r(a+l)cT a ) = E[E p ] J(E, 24,^10 + m T K(l)s£ - Tw (4.19) 

This alternative procedure will increase the calculated value of sur- 
face runoff, and decrease the amount of moisture calculated as infiltra- 
tion. The effect of this difference on the CDF of annual yield will 
be discussed in Chapter 6. 

4.4 Vegetal Equilibrium Hypothesis 

From examining the role of vegetation canopy density in the 
average annual water balance, Eagleson (1978f) observed that for a given 
set of climate and soil parameters and for a given k , Equation (4.3) 

V 

defines s q as a function of M. This relationship is illustrated in 

Figure 4.5 by using four sets of representative soil properties, listed 

in Table 4.1, and the conditions P A = nu. and k =1, for the climates 

A P. v 

A 

of Clinton, Mass, aid Santa Paula, Calif. It can be seen that there 
exists a particular value of M = M for each climate- soil combination at 






which s is a maximum. This point of maximum s corresponds to 
o o 

maximum surface and groundwater runoff, which means, for fixed pre- 
cipitation, that there is minimum evapotranspiration from soil moisture. 
Thu*:, at M * M q , it is expected that the term representing evapo- 
transpiration from soil moisture 

E[E tfCE.M.k^ho) - E[E ] (4.20) 

P A r A 

will be a minimum for a given climate^soil combination. This minimiza- 
tion is seen in Figure 4.6 for the same information as that used in 
Figure 4.5. Note that in Santa Paula, the clay and clay loam soils 
cannot absorb enough water to produce canopy densities greater than 0.4 
and 0.8, respectively, as long as = 1. 

The numerical value of is a matter of some controversy. 
Linacre, et al. (1970), report values of for water plants which range 
from .6 to 2.5 depending upon species. Slatyer (1967, p. 53) states 
that the value of k^ can be greater than one since total evapotranspira- 
tion from a plant community, per unit land area, may exceed that from a 
similar area of bare wet soil due to the larger actual evaporating 
surface area. Kramer (1969, p. 338) however, states that evaporation from 
a plant community never exceeds that from a similar area of wet soil. For 
the present, k^ will not be allowed to exceed one. 

From observations of the relationships presented above, 

Eagleson (1978f) develops the vegetal equilibrium hypothesis mentioned 
in Chapter 2. In the light of the above arguments, this hypothesis says 
that natural vegetal systems of given species will develop a canopy density 
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which produces minimum stress under local climatic conditions. A 
necessary condition for minimum stress is that soil moisture take on 
the largest possible value. Thus, by using this hypothesis, the given 
climate, soil, and plant coefficient determine the equilibrium canopy 
density, M « M q , through the water balance equation, where the soil 
moisture is maximum or, equivalently, where the soil moisture evapo- 
transpiration is a minimum. Figure 4.7 illustrates the relationship 
between the dimensionless evapotranspiration parameter, E, and the 
dimensionless evapotranspiration function, J(E, kv) , for the equilibrium 
condition, M = M q . This plot is obtained by minimizing evapotranspiration 
from soil moisture for a given k^, (k^ = 1 in this case), and E using 
Eagleson’s constant soil column cross-section assumption. The expression 
for E is 


in which 


2$n K(l) TCI) <*> 

TTme 2 
P 


e d+2 
s o 


(4.21) 


P 

n 

♦ ( 1 ) 
m 
d 

*e 


reciprocal of mean time between storms 
porosity 

saturated soil matrix potential 
soil pore size distribution index 
soil diffusivity index 
dimensionless desorption diffusivity 


The other terms have been previously defined. 

Also shown in Figure 4.7 is the M q vs.E relationship for the 

equilibrium condition, P. * nu, . As P* varies from mp , E and thus, 

A " A 
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J(E, M q , k^) will change accordingly, while to the first approximation, 

M will remain constant at M . 

o 

Eagleson (1978f) performs an asymptotic behavior analysis of 
the evapotranspiration function to gain insight into the meaning of the 
parameter, E. The evapotranspiration asymptotes shown in Figure 4.7 
are thereby determined. The intersection of these two asymptotes occurs 
at E = 2 /tt, which separates soil controlled from climate-controlled 
evapotranspiration (Eagleson, 1978d). Thus, low values of E correspond 
to relatively dry, warm climates, while larger values indicate humid 
climates. As can be seen from Figure 4.7, low values of occur for 
low E values, and vegetation densities approaching 1 correspond to a 
large E. 

It can now be seen that observations of canopy cover will 
provide a key to determining the effective properties of a soil for a 
given climate. By using the vegetal equilibrium hypothesis in reverse, 
observations of M q may be used in the water balance to obtain information 
about the soil if the climatic variables are known. 

Figure 4.8, which is a plot of J vs. M q , can be obtained 
directly from the information in Figure 4.7. Thus, from observations of 
vegetation density, the evapotranspiration efficiency, J, can be deter- 
mined. To assure the generality of this relationship, the sensitivity of 
J to its independent parameters is studied. From the expression obtained 
by Eagleson (1978d), the primary parameters other than E and M are: 

lc^ = plant coefficient 

Ah Q accounts for storms which do not fill retention capacity 

8h /e measures effect of surface retention on exfiltration 
o p 

X = parameter of Gamma distribution of storm depth 
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Figure 4.8 

EVAPOTRANSPIRATION FUNCTION FOR NATURAL SYSTEMS 

(lc =1, h = 0.0, K = .5, w/e « 1) 
v o P 



. * 

A 



The plots obtained by varying these parameters over their reported 

ranges are presented in Figures 4,9, 4.10 and 4.11. Figure 4.9, which 

holds k and constant at .5 and 1, respectively, illustrates the 

insensitivity of the evapotranspiration function to changes in Ah Q as 

compared to $h o /e p . By holding 3h Q /e p equal to .1 and k v equal to 

1, k is varied in Figure 4.10. As can be seen, the function changes 

infinitesimally with changes in k. In Figure 4.11, the two variables 

k and Ah are held constant at median values, and 3h / 5 and k are 
o o p v 

allowed to vary. From this analysis, the evapotranspiration function 

is shown to be most sensitive to the two parameters, f3h /e and k . 

op v 

Also shown in Figure 4.11 as dashed lines are the curves obtained 

using the alternate formulation of evapotranspiration. Equation (4.14), 

developed in Section 4.2. In review, this expression was developed 

by accounting for the effect of the vegetated fraction of the soil 

column surface on the vertical flux of the exfiltrating soil moisture 

in Equation (4.6). Expanding the Philip exfiltration equation, which 

was developed for the one-dimensional case of a constant cross-section, 

to two dimensions introduced an inconsistency with the results Philip 

£ 

obtained as explained in Section 4.2. By multiplying the term f g in 

Equation (4.6) by (1-M) , and not adjusting the terms on the right hand 

side of the expression, an infinite exfiltration capacity is obtained 

for the case when M = 1. Although the term (1-M) appears in the 

% 

denominator of several components of the equation (4.14) for bare soil 

storm exfiltration volume, an infinite result is not obtained since the 

total volume of bare soil exfiltration, E , is weighted by (1-M) in 

s 

A8 




Figure 4.9 

SENSITIVITY OF EVAPOTRANSPIRATION FUNCTION TO &h Q / AND Ah c 
(k = 1, K = .5, w/e « 1) 



Figure 4.10 

✓ 

SENSITIVITY OF EVAPOTRANSPIRATION FUNCTION TO K 

(gh /e = .1, k = 1, w/e « 1) 
op v p 
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, Figure 4 . 11 

SENSITIVITY OF EVAPOTRAN SP IRAT ION FUNCTION TO Bh /e AND k 

o p v 

(k = .5, Ah = .05 w/e « 1) 
o F 
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Equation (4.9). The net result of this modification is to raise the 
total evapotranspiration for a given value of M since the rate has 
been increased. This is seen in Figure 4.11 where the dashed lines 
are plotted above the corresponding solid lines for the same vegata- 
tion density. The main problem with this approach, as mentioned 
above, is that the terms on the right hand side of Equation (4.6) 
do not identically go to zero as M approaches 1. If the necessary 
corrections were known, the result would be a reduction in the bare 
soil exfiltration capacity for each value of M. This would lower 
the dashed lines of Figure 4.11. The actual function may therefore 
lie somewhere between th*jte two sets of curves. With this in mind, 
these plots will be used in the following chapters to study the 
validity of the vegetal equilibrium hypothesis, and to determine its 
utility in estimating the effective average areal soil properties of 
a natural watershed. 


52 


j^aagaaifc 


Chapter 5 


METHOD OF APPROACH 


5.1 Vegetal Equilibrium Hypothesis 

Verification of the vegetal equilibrium hypothesis presented 

in Chapter 4 is the first objective of this work. This can be accomplished 

through comparisons of actual data (from watersheds representing various 

types of climates) with the hypothesized relationship illustrated in 

Figure 4.11. In review, the hypothesis states that the vegetation 

denisty seeks that value, M , which maximizes soil moisture. This value 

o 

maximizes water yield and thus, for a given climate and soil, minimizes 
evapotranspiration from soil moisture. Minimum evapotranspiration 
can be translated into a value of evaporation efficiency, J. (i.e. the 
ratio of actual to potential evapotranspiration) leading to the 
relationships previously presented in Figure 4.11. 

The average annual water balance is presented by Eagleson 

(1978a) as 

e[p a ]-e[f Ta j-e[r Sa ]+e[Ec a ]-e[y a ] (5#1) 

which states that average annual precipitation minus average annual 
evapotranspiration will equal the average annual basin yield which is 
composed of surface runoff plus groundwater runoff. When analyzing 
catchment data to calculate average annual actual evapotranspiration, 
mean annual basin yield is subtracted from mean annual precipitation. 
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Although relatively accurate annual precipitation data are readily 
obtained from station records, yield information is only available 
in the form of streamflow records. Therefore, it is necessary that 
the total yield of the catchment studied appear as streamflow. This 
means that the entire groundwater component of yield must be influent 
to the stream channel upstream of the basin mouth. Under such condi- 
tions, most closely approached in humid climates, the total evapo- 
transpiration is equal to precipitation minus streamflow* This 
restriction may lead to overestimating actual evapotranspiration if 
there are losses of yield to ungaged groundwater, or underestimation 
if there is contribution to streamflow of groundwater from adjacent 
watersheds . 

Potential evapotranspiration is estimated by using the 
modified Penman equation (Penman, 1948). The form used here is the 
combination form as presented by Eagleson (1977) 


= q i (l-A) - q b 
P e L e (l"+Y/A) 


in which 

q^ = average rate of insolation 

q^ = average rate of net outgoing long wave radiation 

* 

H = average residual sensible heat flux 
A = shortwave albedo of surface 
P e = mass density of evaporating water 
L £ = latent heat of vaporization 

y/A = atmospheric parameter, a function of atmospheric temperature 
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The above parameters are calculated or estimated as follows: 

qj * <L(4>); from Figure 5.1, where <p m latitude 

A » A (surface structure); from Table 5.1 

1/(1 + Y/A) * f(T^); from Figure 5.2, where T^ "average 
annual temperature 

L * 597 cal/gr 
e 

P e * 1 gr/cm 3 

q b - (1 - .8N) [ .245 -.145 x 10“ 10 T^] 

H - q b /(.25 + 1/(1 - S)) 

The necessary climatic variables are available from U.S. 
Weather Bureau publications. They must be averaged over the rainy 
season which is assumed to be identical with the vegetation growing 
season. 

Equation (5.2) gives the average potential evapotranspiration 
rate. The total potential volume is obtained by multiplying e^ by the 
season length as determined from monthly rainfall records. 

With actual and potential evapotranspiration known, the only 
remaining variables needed for comparison with the hypothesis are the 
vegetation species (to obtain k v ) and the canopy density. The canopy 
density is estimated either from aerial photographs, from personal 
observation, or from literature available for the catchment studied. 

In this work, no photographs were available, and it was possible to 
estimate only ranges of density from the information in the literature, 
depending upon each author's method of measurement and interpretation. 

As a result of this uncertainty regarding the type and canopy density of 
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Table 5-1 


Albedo of Natural Surfaces 
(from Ref. [17]) 


Surface 

Albedo, A 

Surface 

Albedo, A 

Water 

0.03-0.40 

Spring wheat 

0.10-0.25 

Black, dry soil 

0.14 

Winter wheat 

0.16-0.23 

Black, moist soil 

0.08 

Winter rye 

0.18-0.23 

Gray, dry soil 

0.25-0.30 

High, dense grass 

0.18-0.20 

Gray, moist soil 

0.10-0.12 

Green grass 

0.26 

Blue, dry loam 

0.23 

Grass dried in sun 

0.19. 

Blue, moist loam 

0.16 

Tops of oak 

00 

V 

o 

Desert loam 

0.29-0.31 

Tops of pine 

0.14 

Yellow sand 

0.35 

Tops of fir 

0.10 

White sand 

0.34-0.40 

Cotton 

0.20-0.22 

River sand 

0.43 

Rice field 

0.12 

Bright, fine sand 

0.37 

Lettuce 

0.22 

Rock 

0.12-0.15 

Beets 

0.18 

Densely urbanized 


Potatoes 

0.19 

areas 

0.15-0.25 

Heather 

0.10 

Snow 

0.40-0.85 



Sea ice 

0.36-0.50 




• aafreflfrimiMfc • 



LONG-TERM AVERAGE POTENTIAL EVAPOTRANSPIRATION 

(from Ref. (15)) 



the vegetation, the other variable which is a function of the surface 
structure, albedo, is subject to error as well. Therefore, the 
catchments studied will be plotted on Figure 4.11 in terms of the 
expected range of both J and M q . 

5.2 Estimation of Effective Soil Properties 

The second goal of this work is to use the hypothesized 
relationship between vegetation canopy density and evapotranspiration 
to estimate effective average areal properties of the soil. 

Three types of parameters are considered: climate, soil 

and vegetation. The climatic and vegetal properties are easily obtained 
from observations; this leaves the four soil parameters, s q , k(l) , n, 
and c to be determined from the derived relationships between climate, 
soil and vegetation. 

The range of values of the porosity, n, is known to be quite 
small, from .25 to about .45, and does not have a large effect on solu- 
tions of the water balance equation. Assuming a value for n leaves the 
soil moisture, intrinsic permeability, and pore disconnectedness index 
as unknowns. To solve for these variables, three equations or relation- 
ships are needed which incorporate the vegetation and climate as well. 
The first relationship is the water balance. Equation (4.3), which 
expresses the soil moisture, s q , as an implicit function of the climate, 
vegetation and soil. The vegetal equilibrium hypothesis provides the 
second relationship between the same three parameters. The third 
expression used is a rather weakly correlated regression between the 
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intrinsic permeability of the soil, k(l), and its pore size distri- 
bution index, ra, presented in Figure 5.3. This expression is 
(Eagleson, Personal Communication) 

k<1) " ( 5lO )2 ‘ 75 * C “ 2 (5 * 3) 

where 

m - 2/ (c-3) (5.4) 

The coefficient of determination of this regression is small 
due to the extreme variability of these parameters in nature. The effect 
of this regression equation on the derived CDF of annual yield will be 
observed in Chapter 6. 

In order to explain the procedure followed in the estimation 
of soil properties, it is necessary to present mathematically' the water 
balance and the vegetal equilibrium hypothesis. The mean annual water 
balance. Equation (4.3), is again 


mp (1 - e ” G_2a r (a+l)cT a ) = E[E ] J(E,M,k v ,h o ) - E[E r J +m T K(l)s£ - Tw 
A P A A 


for 


E[E J/mp < e" G " 2a T(0 + l)o"° 
A A 

Otherwise, 


(5.5a) 


hl. = E[E ] J(E, M, k , h ) + m K(l) s C - Tw (5.5b) 

T A P A vox o 

If the interpretation of surface runoff developed in Section 


4.3 is used, the above equation becomes 


nu (1 - e“ G " 20 r(a + 1 )a~°) = E[E n ] J(E,M,k ,h ) + m K(l)s 0 - Tw 

r A Pa V O T O 

(5.6) 

Although the components and symbols have been defined earlier, 
their full expressions have not all been stated. In the above equations 


G = aK(l) [ (1 + s q ) C /2 - w/K(l) ] 




nn 2 K(l) ¥(1)(1 - s o ) 2 4>,(d, S > 1/3 


6 it 5 in 




(5.7) 


(5.8) 


E[E ] = m m [1 - M(1 - k )] e 
P A V fc b v P 

K(l) = k(l) y w /y w 


where 


y(d 


in which 


or “5 1/2 

— I — 2 — H 

^w ^K(l) 


a = reciprocal of mean storm intensity 
ri = reciprocal of mean storm depth 
d = 2 + 1/m 

6 = reciprocal of mean storm duration 

m = mean number of storms 
v 

m = mean time between storms 


(5.9) 


(5.10) 


(5.11) 


w 


w 


w 

(j>(m) 


= specific weight of water 
= surface tension of water 
= viscosity of water 


= pore shape parameter = 10 


0.66 + 0.55/m 4* 0.14/m 


2 
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The vegetal equilibrium hypothesis states that 


3E[E* ] 

— - =0 at M = M 

9M o 


where, as mentioned before. 


E[E_ ] « evapotranspiration from soil moisture; 

L A 


which is 


E[E t ] = J(E, M, k y , h Q ) E[E ] - E[E r ] 


with 


E[ V ■ ' ” v 


(1 - M)Sl - e 


f- 


-$h/e r[K, Ah_] 

o p o 


roc) 


r pn /e -| 
_ 1 + — °- . P 

1 + Ah 


3h / e -| -K y[K, (Ah + 3h /e )] 


o P 


T(k) 


+ k 


-3h /e r[K, Ak h 1 


f, " "o' ~p ' "v-'o- 

- M l 1 - e — rW 


r 3h /e -i-k y[ic, (Ak h„ + 3h /e )]• 
h . o P v o o P 

"I 1 Ak h J 

u V o J 


r(K) 


Therefore, 


9E[E* ] 
A 
3M 


3E[Ep ] 3J(E, M, k , It 

J(E, M, k h ) — jg- + ECE ] V ° 

A 


9M 


9E[E ] 

!a_ 

9M 


(5.12) 


(5.13) 


(5.14) 


) 


(5.15) 
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For ■ 1, all M sensitivity comes from J(E, M, k v , h Q ), which will then 
have a minimum at M^, in which special case, according to Equations 
(5.11) and (5.12), 

8J(E, M, k , h ) 

9m - V - -- ° • = 0 at M = M q (5.16) 

As discussed in Section 4.4, evapotranspiration efficiency, J, 

can be determined, for a given climate, from observations of vegetation 

density and species by using the vegetal equilibrium hypothesis, Eq. 

(5, 14). The actual procedure for doing this is to pick a value for the 

evaporation parameter, E, and calculate J for different values of M until 

evapotranspiration from soil moisture, Eq. (5.13), is minimized. If 

the vegetation density obtained which minimizes E is not equal to the 

X A 

observed value, E is incremented and a new M is found. For a fixed 
climate, variations in E correspond to variations in the soil properties 
k(l), c, n, and s q . Therefore, what is actually done is seeking the soil 
which produces the observed vegetation canopy density for a specific 
climate. Once this value of evapotranspiration is found, the value of 
E is also known, which is a function only of the soil parameters for a 
given climate. 

With this information in mind, the following procedure is used 
to estimate the average areal effective soil parameters for a given set 
of climatic and vegetal parameters: 

1. A value for n is assumed and k is set * 1 
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2. The above procedure is followed to determine E 

3. The lowest possible value for c, approximately 3.1, is 
picked as an initial value 

4. k(l) is calculated from Equation (5.3) 

5. With these values for the three soil parameters, n, k(l), 
and c, it can be seen from Eq. (4.21) that s q remains as 
the only variable needed for determining E. With E known 
from step 2, s q is calculated 

6. Annual precipitation is calculated via Equations (5.5) 
through (5.10) 

7a. If the annual precipitation from Step 6 is not equal to the 

actual mean rainfall, c is incremented upward from its 

initially low value and Steps 4-6 are repeated 

7b. Due to the approximation introduced by using Equation (5.3), 

the precipitation, P^, calculated in Step 6 may never 

exactly equal the actual mean value, m- , for any value 

A 

of c. P. will approach trip as c is increased, coming to 
A *A 

within AP^ of equality at intermediate c before diverging 
again for large c. For low values of c, the calculated 
k(l) is large, representing a soil with high permeability 
and well connected pores. With evapotranspiration 
specified at the optimum (i.e,, minimum) value, a large 
precipitation is therefore calculated in order to produce 
the inevitably large groundwater yield of the highly porous 
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soil. For large c and small k(l) , the soil is extremely 
impervious and the surface yield will be high. With 
minimum evapotranspiration, a large value for precipitation 
is again needed. Somewhere between these two extremes, a 
set of suitable soil parameters is obtained which gives an 
annual precipitation, P^, which is closest to the actual 
mean, . This relationship is illustrated in Figure 5.4. 
Holding c constant at the value which gives the minimum 
A p A> k(l) is then deviated from regression equation (5.3) 
until another minimum in calculated precipitation is 
reached. If this value is above the mean precipitation, 
c is decreased, if it is below the mean, c is increased. 
Another search is done on k(l) until the minimum precipi- 
tation is found. This step is repeated until the minimum 
calculated precipitation is equal to the mean 
8. If the values obtained for k(l) and c are not consistent 
with the assumed porosity, n is adjusted to a more appro- 
priate value corresponding to a more pervious or impervious 
soil type depending on the values of k(l) and c. Steps 1 
through 7 are repeated. 

The soil parameters obtained from Steps 1-8 are used to 
construct the CDF of annual yield in the same manner as Eagleson (1978g) . 
In this paper, the annual water balance is written as 
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Figure 5.4 

WATER BALANCE SOLUTIONS USING SOIL PROPERTIES FROM EQUATION (5.3) 






(5.17) 


P A- 


+ R 


In order to relate the annual water balance, (5.17), to the 

mean annual water balance, (5.1), Eagleson defines a climatic mean, m £ , 

where * nip and E^, ■ E[E^, ], and expands (5.17) about this point in 

AAA 

a multidimensional Taylor expansion [Hildebrand, 1959, p. 353]. By 
taking expected values of this expansion term by term, neglecting higher 
order terms, and assuming all variances, covariances, and curvatures 
are small, the "first order approximation" of E[Y^] is obtained: 

EIY A 3 = P A “ E T “ R s + R e (5 * 18) 

A A X A S A 8 A 

This allows the use of the mean annual water balance equation to 
calculate annual values by letting the annual precipitation, and thus 
the average annual soil moisture, vary. The CDF of annual yield can 
then be calculated. Comparison of this CDF with that obtained from 
observations of annual streamflow provides the test for the accuracy 
of the estimated parameters, n, k(l) and c. 

Chapter 6 will present the results of this procedure in the 
form of annual CDF's of basin yield, in audition to verification of the 
equilibrium vegetation density hypothesis. 
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Chapter 6 


PRESENTATION OF RESULTS 


6.1 Verification of Vegetal Equilibrium Hypothesis 

The results of the applied methods of analysis explained in 
Chapter 5 are presented in this chapter. The vegetal equilibrium hypo- 
thesis is verified first in order to assure its validity for use in the 
estimation of soil parameters. 

Appendix A presents the individual catchments studied, the 
data used, location of the catchment, the values obtained for potential 
and actual evapotranspiration, vegetation density, and the estimated 
value of J. 

Figure 6.1 presents the agreement of these experimental data 
with the hypothesized theoretical curves of Figure 4.11. As can be 
seen, the dashed curves, which represent the derivation accounting for 
the presence of vegetation at the surface of the soil column in the 
exfiltration equation, provides a better fit for catchments with a 
vegetal canopy density greater than 0.2. This may mean that the presence 
of vegetation has a much greater effect on soil moisture exfiltration 
than previously believed. Although the equation used has serious flaws, 
they may be negligible compared to the possible importance of the 
presence of vegetation. 

Possible reasons for catchments W-4, W-5, and part of W-8 
lying above the curve may be ungaged yield which escapes through 
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.3 




Figure 6.1 

VERIFICATION OF VEGETAL EQUIL 













groundwater aquifers, or flaws in the vegetal equilibrium hypothesis. 
Until it can be determined if all these yields are present in the 
observed streamflows, the vegetal equilibrium hypothesis would seem to 
give a reasonably accurate relationship between vegetation density and 
evapo transpiration. 

Also shown on Figure 6.1 are the results obtained from the 
empirical formulas developed by Czarnowski (1964) and Ritchie and 
Barnett (1971). These functions exhibit the same type of relationship 
between evaporation efficiency and vegetation density, but do not fit 
the observed data quite so well as Equations (4.10) or (4.14). The fact 
that the data for these studies are primarily from agricultural areas, 
where cultivation and irrigation significantly violate the assumption 
of natural watersheds, is a likely reason for the poor observed fit. 

6.2 Estimation of Effective Average Areal Soil Properties 

To determine the accuracy of the procedure described in Section 

5.2, the two catchments studied here will be those studied by Eagleson 

(1978f, g) ; Clinton, Ma. and Santa Paula, Ca. Table 6.1 presents the 

list of necessary input variables (Eagleson, 1978g) and the computer 

program employed is listed in Appendix B. Tables 6.2 and 6.4 list the 

results obtained for the inputs given in Table 6.1. Listed probabilities 

are calculated for given P /m_ using the Poisson model of Eagleson 

A i A 

(1978b). In Clinton, the value of is held constant for the entire 
range of soil moistures, while in Santa Paula, the vegetation density 
is allowed to vary with annual precipitation, as explained by Eagleson 
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(1978g) . These results are obtained using the form of surface runoff 
developed in Section 4.3 and presented in Section 5.2, (Equation 5.6). 
Figures 6.2 and 6.3 present the results in the form of CDF's of annual 
yield. 

Figure 6.2, which represents Santa Paula, also shows the CDF 
obtained by Eagleson (1978g) for his silty-loam soil, which is listed in 
Table 6.3. The soil properties estimated by the algorithm explained in 
Section 5.2 indicated a slightly less permeable soil than the silty loam. 
This soil gives an improved fit over the entire range of CDF values, 
especially in the critical lower tail. 

The results for Clinton are illustrated in Figure 6.3. The 
soil properties obtained in this case indicate, again, a more impermeable 
soil than the silty-loam employed by Eagleson. Although these values 
for k(l) and c. are quite different, the resulting CDF of annual yield is 
indistinguishable from that obtained for siity-loam. To facilitate the 
comparison between the two results. Table 6.5 lists the annual water 
balance components for Clinton, using the silty-loam soil properties. 

Since the eslimated soil properties represent a tighter soil which reduces 
the mobility of moisture, the soil moisture values are higher than for 
the silty-loam. The other major differences between the two soils are 
the values for surface and groundwater runoff. The more permeable silty- 
loam yields a large groundwater component, and a surface runoff component 
which seems unrealistically low for all values of annual precipitation. 

In the case of the estimated soil properties, the surface runoff is 
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Clinton, Mass. 


Santa Paula 

.15 

e.p cm/day 

.273 

3 

m days 
C b 

10.4 

.32 

m days 

r 

1.43 

365 

days 

212 

.5 

K 

.25 

.1 

h cm 
o 

.1 

0 

w/e 

P 

0 

8.4 

T a °c 
A 

13.8 

94 

nip cm 

*a 

54 

1 

k 

V 

1 

.8 

M 

o 

.4 


Table 6.1 

INPUT CLIMATE AND VEGETATION PARAMETERS 
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EFFECTIVE AREA! AVERAGE SOIL PROPERTIES 
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ANNUAL WATER BALANCE COMPONENTS, SANTA PAULA, CA 
ESTIMATED SOIL PROPERTIES. EQUATION (5.6), M = 
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SILTY-LOAM SOIL PROPERTIES. EQUATION (5.6) 


’/. LESS THAN 



FREQUENCY OF ANNUAL BASIN YIELD, SANTA PAULA, CA 



EFFECTIVE AF EAL AVEPAGE SOIL PROPERTIES 
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ANNUAL WATER BALANCE COMPONENTS, CLINTON 
ESTIMATED SOIL PROPERTIES. EQUATION (5.6), 


EFFECTIVE AREAL AVERAGE SOIL PROPERTIES 
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ANNUAL WATER BALANCE COMPONENTS, CLINTON, MA 
SILTY-LOAM SOIL PROPERTIES. EQUATION (5.6), M = 


V. LESS THAN 







greater and the groundwater runoff lower than for the silty-loam soil. 

The identical CDF f s of yield for the two soils can be explained in the 
following manner: 

In this development, the storage of moisture is not taken 
into account, therefore, yield is equal to precipitation minus evapo trans- 
piration. In the Clinton system, evapo transpiration is controlled primar- 
ily by the climate (Eagleson, 1978d), and is relatively insensitive to 
the soil properties except for extreme cases. Thus, for a given precipi- 
tation, evapotranspiration and hence yield, will be the same for different 
types of soil. The only variations occur in the proportioning of yield 
between surface and groundwater runoff. The permeable soil encourages 
gravitational percolation and hence groundwater, while the impermeable 
soil rejects precipitation as surface runoff. 

In Santa Paula, where evapotranspiration is primarily soil 
controlled, the yield is more sensitive to changes in the soil properties, 
and thus there is a difference in the CDF f s for the two different soils. 

In Figure 6.4, the estimated soil properties are used to show 
the effect on the yield CDF of the two methods of handling surface reten- 
tion in calculating surface runoff. As expected, the values obtained for 
yield, using Eq. (5.5) are reduced from those calculated by Eq. (5.6) due 
to the reduction of rainfall excess in favor of surface retention. 

Although the difference between the two equations is not large. Equation 
(5.6) still fits the observed data better in the lower tail. 

Tables 6.6 and 6.7 list the CDF f s obtained for Clinton and 
Santa Paula, using Eq. (5.5). Again, in the case of Clinton, the CDF is 
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ANNUAL WATER BALANCE COMPONENTS, SANTA PAULA, CA 
ESTIMATED SOIL PROPERTIES. EQUATION (5.5), M = 



EFFECTIVE APEAL AVERAGE SOIL PROPEPT 
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ANNUAL WATER BALANCE COMPONENTS, CLINTON 
ESTIMATED SOIL PROPERTIES. EQUATION (5.5), 


identical to that obtained from !'• . (5.6). This can again be attributed 
to the fact that Clinton is primarily climate controlled, and evapotrans- 
piration is held almost constant near the potential regardless of the 
amount of water that is infiltrated or removed as surface runoff. 

In order to study the sensitivity of the results presented 
here to the vegetation density, values of M q that bracket the observed 
values are used in the soil property estimation program. Figure 6.5 
illustrates the results obtained for Santa Paula, which are listed in 
Table 6.8. Inputing an M q of 0.2 generates a set of soil properties 
that produces more yield and less evapotranspiration than the soil obtained 
using an M q of 0.4. By specifying such a low vegetation density, the 
vegetal equilibrium hypothesis used in the water balance produces a low 
value of evaporation efficiency, J (Figure 4.8). This corresponds to an 
annual evapotranspiration considerably below the potential. By reducing 
the evapotranspiration, the yield must be increased for a given precipi- 
tation, as can be seen by Eq. (5.1). 

On the other hand, attempting to input an M q which is greater 

than 0.41 does not give a solution. That is, no soil can be found for 

the Santa Paula climate which will produce a vegetation density much 

larger than the observed value of .4. The climatic variables, e , mp , 

and m , at Santa Paula prohibit the system from sustaining a larger 
b 

vegetation density, and thus a higher evaporation efficiency. If annual 
precipitation is increased, or e p decreased, the resulting increased 
availability of moisture would allow a greater M . 

The same type of results are seen in Figure 6.6 and Table 6.9 
for Clinton. Even though the vegetatic i density is already large, and 
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ANNUAL WATER BALANCE COMPONENTS, SANTA PAULA, CA 
ESTIMATED SOIL PROPERTIES. EQUATION (5.6), M = 
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ESTIMATED SOIL PROPERTIES. EQUATION (5.6) 


EFFECTIVE AFEAL AVERAGE SOIL PROPERTIES 
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ESTIMATED SOIL PROPERTIES. EQUATION (5.6) 


evapotranspiration is near the potential, it is still impossible to find 
a soil which allows an much larger than the observed value of 0.8. 
Again, reduction of M q produces a soil which generates a larger amount 
of annual yield for the same reasons mentioned for Santa Paula. 

On the basis of these comparisons we see the soil properties 
determined from the estimation algorithm describe the behavior of these 
two systems very well through the water balance model. A brief summary, 
and conclusions drawn from these results will be presented in Chapter 7. 

Although the yield CDF's for Clinton derived from varying 
soil properties are identical, the values obtained for the average annual 
soil moisture vary significantly between the silty-loam soil and the 
soil found from the algorithm. Since soil moisture is a state variable, 
it is desirable to be able to verify the accuracy of its prediction. 

One possible method for doing this would be to compare the CDF's of 
surface runoff, rather than total yield. It has been noticed that the 
surface runoff components of the annual water balance are much more sensi- 
tive to changes in soil properties than is the total yield. One problem 
with this, however, is the lack of measurements of surface runoff, al- 
though streamflow in arid climates may actually be composed totally of 
surface runoff. 
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Figure 6.6 

SENSITIVITY OF ANNUAL BASIN YIELD TO VEGETAL CANOPY DENSITY, CLINTON, MA 





Chapter 7 


SUMMARY AND CONCLUSIONS 


A one-dimensional water balance model (Gagleson, 1978a, b,c,d, 
e,f) is employed to parameterize the climate-soil-vegetation relationship 
at the land-air interface. A vegetal equilibrium hypothesis proposed by 
Eagleson (1978f) provides a second relationship between the climate, soil 
and vegetation. 

Improvements are made in the method of calculating the bare 
soil component of evaporation, and in the way surface retention is 
handled. 

The vegetal equilibrium hypothesis is developed, and its use 
in the water balance is explained. The sensitivity of this hypothesis 
to various parameters of the evapotranspiration function is explored. 

It is found that the two parameters to which the system is most sensitive 
are 8h Q /e^, which can be readily evaluated, and k^, whose value is uncer- 
tain. It is believed that is usually equal to one, except in very 
dry climates, where the plants transpire at a rate less than an equiva- 
lent area of bare wet soil. In this work, k is held at its nominal 

v 

value , 1 . 

Reasonable verification of the vegetal equilibrium hypothesis 
is obtained through comparisons of the theoretical relationship between 
density of canopy cover and the evapotranspiration efficiency to data 
obtained from observations in watersheds representing various types of 
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climates. 

An algorithm is derived which searches for the soil properties 
that produce, in a given climate, the level of evapotranspiration deter- 
mined through observations of vegetation density. By using the vegetal 
equilibrium hypothesis, the water balance, and a regression equation relat- 
ing the soil’s intrinsic permeability and pore size distribution index, 
a consistent set of soil properties is found which generates the implied 
evapotranspiration and also satisfies the mean annual water balance. 

This estimation of soil properties produces results, through 
the water balance, in the form of CDF's of annual basin yield, that 
describe the observed behavior of the Clinton and Santa Paula systems 
very well. In both Clinton and Santa Paula, the soils determined were 
slightly less permeable than the silty-loam which Gagleson (1978g) used 
as his best-fitting soil. These soils also produce a more realistic 
(although unverified) surface runoff component than those used by 
Eagleson. 

A remaining important question is the sensitivity of the 
water balance model to the vegetation parameters, M q and k v » Inclusion 
of this analysis was beyond the scope of this study, and it is left as 
an important subject of future work. 

From this summary, the following conclusions may be drawn: 

1. The vegetal equilibrium hypothesis is sufficiently valid 
to justify its use as a supplementary water balance rela- 
tionship between the soil, climate, and the vegetation. 


92 


2. The algorithm for estimating the effective areal soil 
properties works well, producing CDF's of annual yield 
which fit the observed CDF's closely. 

3. It is more accurate to subtract surface retention from 
the volume of infiltrated precipitation at the beginning 
of the rainfall period than from the rainfall excess. 

4. Use of the vegetal equilibrium hypothesis and the soil 
estimation algorithm should facilitate the incorporation 
of the areal variability of soil properties and soil 
moisture dynamics into global climate models. 


4 
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Chapter 6 


RECOMMENDATIONS FOR FUTURE WORK 

Questions remaining and subjects for future study are: 

1. Evaluation of the Philip exfiltration equation for a 
varying soil column cross-section. 

2. Sensitivity of the water balance to vegetation through 
the parameters, and k^. 

3. Development of a procedure for determining the accuracy 
of predicted values of average annual soil moisture. 
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Appendix A 


DATA FOR CATCHMENTS STUDIED IK VERIFICATION OF 
VEGETAL EQUILIBRIUM HYPOTHESIS 
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W-l 


Location: 
Latitude: 
Rainfall: 
Strearaflow; 
Season Length 


Albuquerque, New Mexico 
<j> - 35°N 
P - 4.37 in. 

Q * .2 in. 

* 4 mos. , July - Oct. 


Cloud Cover: N * .37 

Humidity: S * 39.97% 

Temperature: T • 69.61°F 


E„ 


4.37 - .2 


4.17 in. 


Vegetation Density: ** .12 to ,15 


Albedo: A 

e 

P 

J 


.25 to .3 

15.47 in/season to 14.19 in/season 
.27 to .294 


Watershed Conditions: Rough broken rangeland. About 85% is bare. Sparse 

vegetation consists of short grasses, shrubs, and a few small juniper 
and pinion trees. 


Comments: The value for M is estimated directly from the percent bare 

ground, and taking into account the crown spread of the trees. 


Source*: Hydrologic Data for Experimental Watersheds in the United States, 

1967. U.S.D.A. 

* Indicates reference from which vegetation density values are obtained, 
and in some cases, precipitation and streamflow data as well. All other 
data is obtained from U.S. Weather Bureau publications and U.S.G.S. 
reports of surface water resources. 
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W-2 


Location: Cornfield Wash, New Mexico 

Latitude: <p - 35 °N 


Rainfall: P 

Streamflow: Q 

Season Length 
Cloud Cover : N 

Humidity: S 

Temperature: T 


6.29 in. 

function of M 
o 

4 mos., Jv / - Oct. 
.37 

39.97% 

69.61°F 


Vegetation Density: 

M 

o 

- .16 

Q 

« 1.07 in. 

E-, ■ 5.22 in. 

A 


M 

o 

= .24 

Q 

•» .28 in. 

E_ ■ 6.01 in. 
A 

Albedo: 

A 

= 

m 

CM 

a 

to 

.30 


e 

P 

- 

15.47 

in/ season to 

14.19 in/season 

M * .16 
o 

J 

as 

.34 

to 

.37 

M = .24 
o 

J 

= 

.39 

to 

.42 


Watershed Conditions: The dominant vegetation is galleta grass. Remain- 

ing areas have a mixture of other grasses, Russian thistle, and big sage- 
brush in small upland drainages. 

Comments: Runoff data was recorded as a function of percent bare soil in 

the paper used as the source, therefore, the calculation of E T gives two 

A 

values, one for each and Q data pair. Vegetation density values were 
recorded for each Vi - lue of percent bare soil, and the two extreme values 
were used here. 


100 


W-2 (continued) 


Source: 


F. A. Branson and J. B. Owen, "Plant Cover, Runoff and Sediment 
Yield Relationships on Mancos Shale in Western Colorado," W.R.R. , 
6 ( 3 ), 1979 . 
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W-3, W-4, W-5 


Location: 

Tombstone, Ariz. 




Latitude: 

4 - 32' 

®N 




Rainfall: 

P W-3 

8.65 in 

*r 4 

- 8.01 

in. 

Streamflow: 

%-3 ~ 

.64 in. 

V-3 




P W-4 " 

8.65 in. 

e t a 

- 8.44 

in. 


Q W-4 “ 

.21 in. 

\ f-4 



, 

P W-5 

8.65 in. 

e t a 

- 7.56 

in. 


Q W-5 * 

1.09 in. 

\-5 



Season Length 

m 

3 mos., July - Sept. 




Cloud Cover: 

N - 

.35 




Humidity : 

S - 

.4467 




Temperature: 

T - 

82.17°F 





Vegetation Density: 


M 

.35 

to 

.4 

° W-3 




M 

.25 

to 

.3 

°W-4 




M * 

.2 

to 

.25 

°W-5 





Albedo: 


A « 

.24 

to 

.30 

e * 
P 

13.45 in/season 

to 

12.14 

J W-3 “ 

.60 

to 

.66 

J W-4 " 

.63 

to 

.70 

J W-5 * 

.56 

to 

.62 

All watersheds 

have cover of desert 

shrubs 


(whitehom, creosote bush, tarbush) with an understory of grass (black 
grama, tobosa grass, blue grama, sideoats grama, and curly mesquite grass). 
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W-3: Entire area covered by shrubs with 38% crown spread. M = .35 to .40. 

o 

W-4: 78% of area covered by shrubs with crown spread of 30%. Remaining 

22% covered with grass with .2% basal area. M s .25 to .3. 

o 

W-5: Shrub canopy approximately 20%. Remaining area covered by grass 

with .2% basal area. M * .2 to .25. 

o 


Comments: The three watersheds are all sub-catchments of a larger catch~ 
ment. Therefore, while vegetation densities and streamflow vary slightly, 
the annual climatic properties are all the same. 


W-6 


Location: 

Flagstaff, Arizona 



Latitude 

* = 35°N 



Rainfall: 

P » 12.38 in. 

E - 11.98 
A 

in. 

Streamflow: 

Q «* .4 in. 



Season Length 

* 7 mos., July Jan. 



Cloud Cover: 

N = .4 


. 

Humidity: 

S * .52 



Temperature: 

T = 47.13°F 



Vegetation Density: M q = .3 to .35 



Albedo : 

A = .2 

to 

.25 


e ** 20.63 in/ season to 

P 

18.95 in/ season 


J = .58 

to 

.63 


Watershed Conditions: The terrain is undulating uplands dissected by many 

small drainages. The vegetation is mainly upper pinion juniper woodland 
with a sparse understory of grasses. 

Source: Brown, H.W. , "Characteristics of Recession Flows from Small Water- 

sheds in a Semiarid Region," W.R.R. , 1(4), 1965. 


104 






Location: 

Badger Wash, 

Colorado 

Latitude: 

<j> « 38°N 


Rainfall: 

P = 4.69 in. 


Streamflow: 

Q = function 

of M 
o 

Season Length 

= 6 mos . , 

August. - 

Cloud Cover: 

N = .5 


Humidity : 

S = .4817 


Temperature: 

T « 47.8°F 


Vegetation Density: M q = . 

13, 


M = . 
o 

26, 

Albedo : 

A = 



e = 16 

P 

M = .13: 
o 

M = . 26 : 
o 


f-7 


Jan. 


Q = .96 in.. 


E t - 3.73 
A 

in 

Q = .35 in. , 


E - 4.34 

A 

in 

25 

to 

.30 


04 in/season 

to 

14.66 in/season 

J = .23 

to 

.25 


J = .27 

to 

.30 



Watershed Conditions: The catchment is in a semiarid area with pre- 

dominantly desert-type shrubs. 

Comments: This data was obtained in the same way as that for W-2. Thus, 

the values for J are presented in the same way . 

This watershed is located in an area where there is considerable 
snowfall. The model used in this work does not account for snowmelt in any 
way, and only works with yield resulting from precipitation in the form of 
rainfall. Therefore, if the yield measurement includes runoff from snow- 
melt, the value of precipitation used here is not large enough to account 
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W-7 (continued) 

for that much strearoflow, and the resulting calculated value of actual 
evapotranspiration is too small. It would not be surprising then if the 
value plotted for J vs. M q is below the hypothesized curve. 

Source: Branson, F. A. and J. B. Owen, "Plant Cover, Runoff, and Sediment 

Yield Relationships on Mancos Shale in Western Colorado," W.R.R. , 6(3), 
1979. 


Location: 

Santa Paula, California 


Latitude: 

- 

34.4°N 


Rainfall: 

P = 

21.26 in. 

E * 14.41 in. 

t a 

Streamflow: 

Q - 

6.85 in. 


Season Length: 

m 

7 mos., Oct. - Apr. 


Cloid Cover: 

N * 

.37 


Humidity: 

s = 

.6897 


Temperature: 

T - 

53.06°F 


Vegetation Density: 

M = .35 

o 

to .5 

Albedo: 


A = .2 

to .32 


e = 21.23 In/season to 16.73 in/season 

P 

J = .68 to .86 

Watershed Conditions: Fairly rugged terrain with wide variation of vege- 

tation type. Dominant species are desert- type shrubs which are common in 
Southern California mountain ranges. 

Source: 1) Eagleson, P. S. , "Climate, Soil and Vegetation," Parts 1-7, 

W.R.R., 14(5), Oct. 1978. 


2) On-site observations. 


W-9, W-10 


Location: 

Latitude: 


Chickasha, Oklahoma 
<J> - 35°N 


Rainfall: 

Streamflow: 

P„ . * 23.52 in. 
W-9 

Q w _ 9 “ 1*12 in. 


V,' 

22.40 in. 


P tI » 23.52 in. 

W-1U 

^w-io “ 3,77 in ‘ 


\-10 

- 19.75 in. 

Season Length: 

■ 7 mos . , Apr . - 

Oct. 



Cloud Cover: 

N ■ .47 




Humidity : 

S - 67% 




Temperature: 

T * 70. 61°F 




Vegetation Density: M = 

°W-9 

.45 

to 

.57 


M 

.2 

to 

.3 


°W-10 




Albedo : 

A 

.18 

to 

.24 


e = 

P 

28.80 in/ 

season to 

26.09 in/season 

W-9: 

J 

00 

r^. 

• 

to 

.86 

W-10: 

J 

.69 

to 

.76 


Watershed Conditions: The vegetation of both catchments consists of native 

grasses (buffalo grass, blue grama, little bluestem). Values of M q are 
interpreted from radiation shielding values obtained from average values of 

leaf area index and percent mulch cover. The equation used is (2) 

R -,4(L. .+2.5M) 

M = 1 - = e Al 

0 R no 
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♦wJM| 



where 

R = net radiation reaching the soil surface 
ns 

R no = n t .. radiation above plant canopy 

L. . = leaf area index 
Ai 

M = fraction of surface covered by Mulch 

Comments: It is reported in the source paper that W-10 is constantly 

overgrazed, thus, it is likely that the value obtained for M q is unnaturally 
small, and the plotted position of this catchment will be above the hypo- 
thesized curve. 

Source: 1) Hydrologic Data for Experimental Agricultural Watersheds in the 

U.S. 1976. U.S.D.A. 

2) J. T. Ritchie, E. D. Rhoades and C. W. Richardson, "Calculating 
Evaporation from Native Grassland Watersheds," Transactions of 
the A. S.C.E . , Aug. 1976. 



W-ll 


Location: Clinton, Massachusetts 

Latitude: 4> * 42.50°N 


Rainfall: 
Streamflow: 
Season Length: 
Cloud Cover: 
Humidity : 
Temperature: 


P * 4J.82 in. 
Q = 21.81 in. 

= 12 mos. 

N - .35 
S = .70 
T = 47.12°F 


Vegetation Density: M q = 

Albedo : A = 


.8 

.25 


E = 22.01 in. 

*■ k 


to . 9 

to . 30 


e = 24.25 in/ season to 21.64 in/season 

P 

J = .91 to 1.02 


Watershed Conditions: No specific conditions are available, only the range 

of vegetation density. 


Source: 1} Eagleson, P. S., "Climate, Soil and Vegetation," Parts 1-7, 

W.R.R. , 14(5), Oct. 1975. 

2) Visual observations of nearby watersheds. 
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Appendix B 


FORTRAN PROGRAM FOR ESTIMATION OF SOIL PROPERTIES 


o THIS PROGRAM CALCULATES EFFECTIVE AREAL AVERAGE SOIL PROPERTIES. WHEN 
c THE SOIL PROPERTIES ARE VARIED USING A REGRESSION EQUATION , CALCULATED 
c PRECIPITATION, Pa, REACHES A MINIMUM AT A MEDIAN VALUE OF c. THE PARA- 
c METER, k(1), IS THEN DEVIATED FROM THE REGRESSION UNTIL ANOTHER MINI- 
c MUM Pa IS FOUND. DEPENDING ON WHETHER THIS VALUE FOR Pa IS ABOVE OR 
c BELOW THE KNOWN VALUE OF mpa, THE PARAMETER c, IS INCREMENTED UP OR 
c DOWN, AND k( 1 ) IS SEARCHED AGAIN UNTIL ANOTHER MINIMUM IS REACHED, 
c THIS INCREMENTATION AND SEARCHING IS CONTINUED UNTIL THE MIMIMUN Pa 
c FOUND IS EQUAL TO mpa. 


integer change , f tm , cfbl , runs, number , mo n, iter 

real *8 mnu,p1 

real mtb,mtr,mh,mpa 

real mi,mo,m,n,nu,k1 ,k2 

c DIMENSIONLESS INFILTRATION DIFFUSIVITY 

fii(d,so)= 1 ./(d»(1 .-so)»»(1 .45-.0375*d)+5./3.) 

c PORE SHAPE PARAMETER 

fi(em)=10.**( ,66+.55/em+.14/em**2.) 

print,'Input parameters in the correct units.' 
print, 'ep, cm/d ay mtb,days mtr,days tau,days kappa,-.' 
print, 'ho, cm w/ep,- ta, degrees C. ' 
input , epr , mtb , mtr , tau , ak, ho , wep , ta 

pistol=1 

c IF pistol=1 , THE ARRAY OF FACTORIALS IN THE CDF SUBROUTINE HAS NOT 
c BEEN CALCULATED YET. ONCE pistol=2, THE FACTORIALS HAVE BEEN STORED 
c AND THE LINES WHICH DO THIS CALCULATION ARE THEN SKIPPED. 

5 print, 'Input mpa, cm kv,- Mo,- n,- .' 

input , mpa , akv , mo , n 

print, 'For annually varying Mo, type 1, for constant Mo, type 2' 

input, mon 

if ( mon . eq . 0 ) st op 
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changes 1 

o IF changes 1, SOIL PROPERTIES ARE NOT YET DETERMINED. IF change=2, THE 
c SOIL PROPERTIES HAVE BEEN DETERMINED AND ONLY THOSE STEPS NEEDED FOR 
c DERIVING THE CDF OF THE WATER BALANCE COMPONENTS ARE USED. 


runssl 

c IF runssl, THIS IS THE FIRST SET OF SOIL PROPERTIES USED, AND NO COM- 
c PARISON OF CALCULATED Pa IS POSSIBLE. IF runas2, THE NEW VALUE OF Pa 
c IS COMPARED TO THE OLD VALUE TO SEE IF A MINIMUM HAS BEEN REACHED. 

cfblsl 

c IF cfblsl, THIS IS THE FIRST DEVIATION OF k(1) FROM THE REGRESSION 
c AND NO COMPARISON OF Pa IS DONE. IF cfbls2, THE Pa CALULATED WITH 
c THIS k(1) IS COMPARED TO THAT CALCULATED USING THE PREVIOUS k(1) TO 
c SEE IF THE SECOND MINIMUN HAS BEEN REACHED. 


ftm=1 

c IF ftmsl, THE SECOND MIMIMUM HAS JUST BEEN FOUND, BUT IF THIS MINIMUM 
c Pajfapa , c MUST BE CHANGED AND THE ENTIRE PROCESS MUST BE REPEATED, 
c THE VALUE OF THE DIFFERENCE BETWEEN THE MINIMUM Pa AND mpa, awbal, IS 
c PRESERVED AND COMPARED TO THE NEXT ONE OBTAINED. THIS COMPARISON IS 
o SIGNALED WHEN ftm=2. WHEN awbal < .001, THE SOIL PROPERTIES HAVE BEEN 
o FOUND. 

c SET INITIAL VALUES 

p1=0.0 
30 = 0.0 
dcs=.1 
dics=. 1 


c *************************** 

c COMPUTE WATER CONSTANTS 
c sut=SURFACE TENSION 

c nusVISCOSITY 

c gamsw=SPECIFIC WEIGHT 


113 


call WATCN(ta,sut,nu f gamsw) 


c COMPUTE CLIMATIC PARAMETERS 

deltasl ,/mtr 
mh=mpa/ ( tau/ ( mt b+zutr ) ) 
mnu=tau/(mtb+mtr) 
mlsmh/mtr 
eta=1 ,/mh 
alphas 1 ./ml 
pi=3. 14159 
betas 1 ,/mtb 

epa=epr*tau*mtb/(mtb+mtr) 
alsak/mh 
alhsal*ho 
bhesbeta*ho/epr 
if(ho.eq.0.0)goto 10 
blesbeta/ ( al *epr ) 
goto 20 
10 blesO.O 

20 alkhsalh*akv 

blkesble/akv 

30 if(change.eq.1)goto 40 


print 

print 

print , • so E Mo Pa/mpa Ya/mpa 

do 400 is 1,45 
sosso+,02 
esecnst*so**d2 
fiidsfii(di,so) 


RSA 


40 


if ( change. eq. 2) goto 45 
goto 50 

45 if(mon.eq.2)goto 60 

goto 55 

c TO SPEED UP THE SEARCH FOR THE VALUE OF e THAT MINIMIZES ETA AT THE 
c OBSERVED Mo, e AND m ARE GIVEN INITIAL VALUES DEPENDING ON THE VAL- 
c UE OF Mo. BY PICKING A VALUE FOR e, THE m THAT MINIMIZES ETA CAN BE 
c FOUND. IF THIS m*Mo, ANOTHER e IS PICKED UNTIL msMo. 

50 if(mo.ge..2)es.3 

if(mo.ge..3)e=.5 
if(mo.ge. .4)e=1 . 
if (mo.ge. .6) es 3. 
if (mo.ge. .7)e=6. 
if(mo.ge..8)e=10. 
if (mo.ge. . 9)e=20 . 

55 if(e.ge..01)bm=.1 

if ( e . ge . • 1 ) bm= . 4 
if(e.ge.1 .)bm=.6 
if(e.ge.10.)bms.9 
if(mo.lt..4)de=.01 
if(mo.ge..4)de=.1 
numbers 1 

60 iter= 1 

dm=.01 

70 bmkvsbm*akv 


C COMPUTE EVAPOTRAN SP IRATION PARAMETERS, B & C. 

bs((l ,-bm)/(1 .+bmkv)+(bmkv*bm)/(2.*(1 .+bmkv)**2.) ) 
if (bmkv.eq.O.O)goto 80 
csl ,/(2. # (bmkv # bmkv)) 
goto 90 
80 csI.elO 

90 be=b*e 
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ce=amin1 ( c*e , 80 . ) 


gamk=gamt ( ak, alh ) /gamma ( ak) 

gamkl =gam t ( ak , alh+bhe )/gamma ( ak) 

gambe=gamt( 1 . 5 , be ) 

gamce=gamt(1 .5,ce) 

gam kv=gamt ( ak , al kh )/ gamma ( ak) 

gamkvl=gam t ( ak , ( alkh+bhe ) ) /gamma ( ak) 


c COMPUTE ANNUAL EVAPORATION FROM f’IRFACE RETENTION 

era=epr/beta*( (1 .-bm)*( 1 ,-exp(-bhe)*(1 .-gamk)-(1 .+ble) , *(-ak) 
& *gamkl) 

4 +bmkv*(1 .-exp(-bhe)*(1 .-gamkv)-(1 ,+ble)**(-ak)*gamkvl))*mnu 

eram=era 


c COMPUTE INTERSTORM BARE SOIL EVAPORATION 

es j=gamk-( 1 .+ble) **(-ak) *gamkl»exp( -be)+ 

4 (1 .-gamk)*(1 ,-exp(-be-bhe) # ( 1 .+bmkv+sqrt(2.»b)*e-wep) 

4 +exp ( -ce-bhe ) * ( bmkv+sqrt ( 2 . *c ) *e ) 

4 +sqrt ( 2 . *e ) *exp ( -bhe ) * ( gam ce-gambe ) ) 

4 +(1 .+ble)**(-ak)*gamkl*(sqrt(2.*e)*(gamce-gambe) 

4 +exp(-ce)*(bmkv+sqrt(2.*c)*e) 

4 -exp ( -be) *( bmkv+sqrt ( 2 . *b) •e-wep ) ) 


c COMPUTE EVAPOTRANSPIRATION FUNCTION 

hj=1 ./( 1-bm+bmkv)*( (1-bm)*esj+bmkv) 
ETN=hj*( 1 ,-bm+bmkv) 


95 


if (change. eq. 2) goto 95 
goto 100 

if (mon.eq.2)goto 160 


c THESE LINES FIND THE M THAT MINIMIZES ETA. 

c IF iter=1 , IT IS THE FIRST TIME THROUGH AND NO COMPARISON IS MADE 
100 if(iter.eq.1)goto 120 

if (abs(dm). It.. 000001 )goto 150 
if (ETN.gt.EIMIN)goto 110 
goto 120 

110 bm=bm-1.5*dm 

dm=-.5*dm 
goto 130 
120 ETMIN=ETN 

bmin=bm 
iter=2 
bm=bm+dm 

130 if (bm) 140 ,70, 70 

140 bmr .1»(bm-dm) 

qrq+1 

if(q.lt.4)goto 70 

c AT THIS POINT, NO Mo CAN BE FOUND THAT IS GREATER THAN 0,AND NEW PAR- 
c AMETERS MUST BE INPUT. 

goto 395 

150 bm=bmin 

ETN=ETMIN 

160 if (change. eq. 2) goto 230 

c THESE LINES FIND THE E CORRESPONDING TO THE GIVEN Mo. 
c IF number= 1 , IT IS THE FIRST TIME THROUGH AND NO COMPARISON IS MADE. 
diff=mo-bm 

if(abs(diff) . It. .0001)goto 200 
if (number. eq. 1 ) goto 170 
if ( dif f »dif fold. le. 0.0) goto 190 
if (number. eq. 2) goto 180 
170 if (diff.lt. 0.0)de=-1 .*de 
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numbers 2 

180 diffoldsdiff 
e=e+de 
goto 60 

190 de=-de*.5 

diffoldsdiff 
ese+de 
goto 60 

200 continue 


AT THIS POINT, THE VALUE OF e HAS BEEN DETERMINED, AND SOIL PROPER- 
c TIES ARE NOW SEARCHED. 

cs=4. 

210 ms2./(cs-3.) 

ficsfi(m) 
dEs2.+1 ,/m 
discs- 1 ,/m-1 . 
d2sdE+2. 
fiedsfie(dE) 

c REGRESSION EQUATION 

k1s(m/51 2.7) **2.75 

k2=k1 

dklskl/10. 

220 continue 

bk1sk1*gamsw/nu 

sil ssqrt ( n/( kl *f ic) ) *sut/gamsw 
sigcsn*eta**2 . *bk1 »si1 / ( pi*m*del ta) *72000 . 
ecnsts2 . *beta*n*bk1 *si1 *fied/( pi*m*epr**2 . ) *86400 . 

c SOIL MOISTURE IS CALCULATED. 

sos(e/ecnst)**( 1 ,/d2) 
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fiid=fii(di,so) 


C COMPUTE WATER BALANCE 

c COMPUTE ANNUAL EV A POTR AN SPIRATION 

230 ETA=ETN*epa 


C COMPUTE ANNUAL GROUNDWATER RUNOFF 
RGA=tau*bk1*so**cs # 86400 

sigrf = ( sigc »f iid«( 1 . -so) »*2 . ) ••. 33333 
g=alpha*bk1 *86400*. 5*0 .+so**cs) 
blop=g+2 . *sigrf 
if(blop.gt.85.)blop=85. 

blip=exp( -blop) *gamraa(sigrf+1 . ) •sigrf *•( -sigrf) 
if (blip.gt..95)blip=.95 


COMPUTE PRECIPITATION, YIELD, RUNOFF 
Pa=(ETA+RGA)/(1 .-blip) 

RSA=blip*Pa 
YasRSA+RGA 

if(change.eq.2)goto 380 
awbal=Pa-mpa 

c NOTE-awbal IS THE DIFFERENCE BETWEEN CALCULATED Pa AND KNOWN mpa. THE 
c FOLLOWING LINES WILL PERFORM THE SEARCH FOR SOIL PROPERTIES THAT PRO- 
c DUCE Pa=mpa . 

if (cfbl.eq.2)goto 260 
if (ftm.eq.2)goto 280 
if(runs.eq.1)goto 250 

c THESE LINES PERFORM THE FIRST MINIMIZATION WHICH ADHERES TO THE RE- 
c GRESSION EQUATION. 
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if(abs(dcs) .lt..001)goto 260 
if(awbal.gt.awbol)goto 240 
goto 250 

240 csscs-1 .5*dcs 

dcs=-.5*dcs 
goto 210 

250 awbolsawbal 

cs=cs+dcs 
runs=2 
goto 210 

c THESE LINES PERFORM THE SECOND MINIMIZATION WHICH HOLDS c CONSTANT 
c AND DEVIATES k(1) FROM THE REGRESSION. 

260 if(cfbl.eq.2)goto 270 

if(abs(awbal).lt..001)goto 320 
if(cfbl.eq.1)goto 280 

270 if(abs(dk1).lt.k2/1000.)goto 320 

if(awbal.gt.awbol)gotc 290 

280 awbolsawbal 

klrkl-dkl 

c SINCE k(1) VARIES BY ORDERS OF MAGNITUDE, dkl MUST BE REDEFINED IF 
c k( 1 ) GETS TOO BIG OR SMALL. 

if(k1.1t.k2/9)goto 300 
if(k1.gt.k2«9)goto 310 
cfbl=2 
goto 220 

290 k1=k1+1.5*dk1 

dk1=-.5*dk1 
goto 220 

300 dk1=dk1/10. 

k2=k1 
goto 220 
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310 dk1rdk1»10 
k2=k1 
goto 220 

c THESE LINES PERFORM THE SEARCH ON c WHICH LOCATES THE MINIMUN Pa FROM 
c THE ABOVE PROCEDURE WHICH EQUALS mpa. 

320 cfbl=1 

if ( abs ( awbal ) . 1 1 . . 0 1 ) goto ) 

if(ftm.eq.1)goto 330 
if(awbal*awbold.lt.0.0)goto 350 
goto 340 

330 if(awbal.gt.0.0)dics=-1.»dics 

ftm=2 

340 awbolds awbal 

CS=CS+diC3 

goto 210 

350 dicss-dics».5 

awbold=awbal 
cs=cs+dica 
goto 210 


360 print, • AVERAGE EFFECTIVE PARAMETERS' 

print 370,e, so, kl ,cs 

370 format( '£*' ,f6.3,2x, 'so*' ,f5.3,2x, 'k(l)** ,el6.7,2x, 'os' ,f6.3) 

change=2 
so =0.0 
goto 30 

380 yl sYa/mpa 
plsPa/mpa 

if(p1.1t..2)goto 395 

c COMPUTE CDF OF PRECIPITATION 

cai 1 PROBZ ( mnu , p 1 , prob , ak ) 
if (prob.lt.. 009 )goto 395 


121 








o o bo 


if(prob.gt..99)goto i»10 
print 390,so,e,bm,p1 ,y1 , RSA , RG A , ETA , orob 
390 formafc(f4.2,3x,f6.3,3x,f5.3,3x,f7.4,3x,f7.4,3x,f7.4,3x,f7.3,3x 

4 ,f6.3,3x,f7.4) 

395 continue 

400 continue 

410 goto 5 

end 


c 


THIS FUNCTION COMPUTES THE INCOMPLETE GAMMA FUNCTION, 
function gamt(a,x) 
if (x.eq.O)goto 40 
if(x.gt.100)goto 50 
sum= 1 ./a 
an=1 .0 
oldssum 

33 old=old # x/(a+an) 

if (old/sum-1 .e-6)20,10,10 

10 an=an+1 . 

sum=sum+old 
if(an-300. >33,33,12 

12 continue 

20 xxx= ( a*alog( x) +alog( sum) -x) 

if (xxx.lt. -80. )goto 40 
gamt=(exp(xxx)) 
goto 60 

40 g amt =0.0 

goto 50 

50 g&AC-gamma(a) 

60 return 

end 


This function computes the gamma function by a Stirling approx. 
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function gamma(y) 
xry+l . 

pi*3. 14159 
stir1rl./(i2.»x) 
stir2=1./(288.»x»»2.) 
stir3=-139./(5l840.»x«»3.) 
stir4=-571 ./ (2488320. »x»»4.) 
atirs1+stir1+stir2+stir3+stir4 
gamnasexp ( -x ) •x** ( x- . 5 ) •sqrt ( 2 . # pi ) *stir/y 
end 


subroutine WATCN(ta,sut,nu,gamsw) 
real nu,nut 

dimension sutt(1 1) ,nut(1 1) ,gamst(1 1) 

data sutt/75. 6, 7^.9,74.2,73. 5, 72. 0,72. 1,71 .4, 70. 7, 70.0, 

& 69.3,68.6/ 

data nut/17 .93e-3 , 15 . l8e-3 , 1 3 .09e-3 , 1 1 .44e-3 , 1 0 .08e-3 ,8 . 94e-3 , 
& 8 . e-3 , 7 .2e-3 , 6 . 53e-3 , 5 . 97e-3 ,5 . 94e-3/ 

data gamst/0. 99987, 0.9999999, 0.99973, 0.99913, 0.99823, 0.99708 
& ,0.99568,0.99406,0.99225,0.99025,0.98807/ 

if(ta.gt.50.)go to 10 
ita*ifix(ta».2)+1 
frac=ta- float ( if ix( ta ) ) 
italsita+1 

sut = ( sut t ( it a 1 ) -sut t ( i t a ) ) *0 . 2*f r ac+ sut t ( i t a ) 

nus ( nut ( ital ) -nut ( ita) ) •0.2*frac+nut ( ita) 

gamsw=( ( gamst ( ital ) -gamst ( ita) ) *.2»frac+gamst ( ita) ) *980 . 

return 

10 sutssutt(ll) 
nusnut(ll) 
gamswsgamst(1 1) 
return 
end 


c DIMENSIONLESS EXFILTRATION DIFFUSIVITY 
function fie(d) 
dimension y(6) 

data y/0. 18, 0. 1 1,0. 077,0. 056, 0.01*4, 0.03V 
if(d.gt.7.) goto 10 
x=d-1 . 

if(x.lt,.1 .)x=1 . 

islfix(x) 

frac=x-float(i) 

yl =alog(y(i)) 

y2salog(y(i+1)) 

fie=exp( ( y2-y1 ) *frac+y1 ) 

return 

10 fies.03*» 

return 
end 


subroutine PR0BZ(mnu,p1 ,prob,ak) 
c THIS PROGRAM COMPUTES THE CDF OF NORMALIZED PRECIPITATION, 
real *8 fac(500) 

real*8 x,a,dlog,gama,gamlid,eps 
real*8 m,k,w,t,z,zl,zu,inz 
real *8 pi ,mnu 

real*8 xold,xsum,sum1 ,sum2 , tot, vtot, void, vnew 
integer v,vm,vmax 

if(pistol.eq.2) goto 301 
do 300 j=1 ,500 
vtot=0.0d0 
do 700 iv=1 , j • 

700 vtotsvtot+dlog(dble(float(iv)) ) 

fac( j)=vtot 

300 continue 

301 continue 
eps=1 .e-5 
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pistol=2 

wsmnu 

t=1. 

ksak 

c INITIALIZING VALUES 
c 

m=w*t 

z=p1 

vm=ifix(sngl(m)) 
vmaxsif ix( sngl(3 . *m) ) 

3 x=m*k*z 

ii=0 

Jj = 1 

suml =O.OdO 
sum2s0.0d0 
13 v=vm-ii 

if(v.eq.O)goto 500 
23 if (v.eq.vmax)goto 600 

c 
c 

c COMPUTE LOG INCOMPLETE GAMMA DISTRUBITION 

a=dble(float(v))*k 
xold=1 .OdO/a 
xsum=1 .OdO/a 
i=1 

100 xold=(xold/(a+i))*x 

xsum=xsum+xold 

if((xold/xsum) .le.eps)gotc 200 
i=i+1 
goto 100 
200 continue 

call algama(a,gamm,ier) 
gamlid=a*dlog(x)-x+dlog(xsum)-dble(ganim) 
c 

c COMPUTE THE SUMMATION OF ALL V TERMS 
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vold=dble(float(v))«dlog(m)-fac(v)+gamlid-m 

if ( void . le . -85 . ) vold=-85 . 

vnewsdexp ( void ) 

if ( v.gt.vm)goto 800 

suml =sum1 +vnew 

if((vnew/sum1).le.eps)goto 500 

ii=ii+1 

goto 13 

500 vsvm+jj 
goto 23 

800 sum2=sum2+vnew 

if((vnew/sum2).le.eps)goto 600 
jj=jj+1 
goto 500 
c 

c COMPUTE CDF OF NOHMALIZED PRECIPITATION 
c 

600 if(m.gt.85.)m=85. 

prob=sum1 +sum2+dexp ( -m) 

return 

end 


